Adapted from SciMLBenchmarks.jl multi-language wrapper benchmark.
# Imports
using LinearAlgebra, Statistics
using OrdinaryDiffEq, DiffEqDevTools, ParameterizedFunctions, Plots
using ODEInterface, ODEInterfaceDiffEq
using Sundials
using SciPyDiffEq
using deSolveDiffEq
using MATLABDiffEq
using LSODA
using ProbNumDiffEq
# Plotting theme
theme(:dao;
linewidth=8,
linealpha=0.7,
markersize=5,
markerstrokewidth=0.5,
legend=:outerright,
)
# Constants used throughout this benchmark
const DENSE = false # used to decide if we smooth or not
const SAVE_EVERYSTEP = false
false
# Utility to benchmark the ProbNumDiffEq.jl solvers
function MyWorkPrecision(prob, abstols, reltols, alg; appxsol, kwargs...)
println(alg)
times, errors = [], []
for (atol, rtol) in zip(abstols, reltols)
sol_call() = solve(
prob, alg; abstol=atol*1e2, reltol=rtol*1e2,
dense=DENSE, save_everystep=SAVE_EVERYSTEP, kwargs...)
sol = sol_call()
tbest = 10000
for i in 1:10
t = @elapsed sol_call()
tbest = min(tbest, t)
end
err = mean(abs.(appxsol(prob.tspan[2]) - sol.u[end]))
push!(times, tbest)
push!(errors, err)
#@info "Info" atol rtol tbest err
end
return errors, times
end
MyWorkPrecision (generic function with 1 method)
# Problem definition and reference solution
lv = @ode_def LotkaVolterra begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a b c d
p = [1.5,1,3,1]
tspan = (0.0,10.0)
u0 = [1.0,1.0]
prob = ODEProblem(lv,u0,tspan,p)
sol = solve(remake(prob, u0=big.(prob.u0), tspan=big.(prob.tspan)),
Vern9(), abstol=1/10^20, reltol=1/10^20)
test_sol = TestSolution(sol)
plot(sol, title="Lotka-Volterra Solution", legend=false)
# Non-probabilistic solvers:
_setups = [
"Julia: Tsit5" => Dict(:alg=>Tsit5())
"Julia: Vern9" => Dict(:alg=>Vern9())
"Julia: RadauIIA5" => Dict(:alg=>RadauIIA5())
"Hairer: dopri5" => Dict(:alg=>ODEInterfaceDiffEq.dopri5())
"Hairer: rodas" => Dict(:alg=>ODEInterfaceDiffEq.rodas())
"MATLAB: ode45" => Dict(:alg=>MATLABDiffEq.ode45())
"MATLAB: ode113" => Dict(:alg=>MATLABDiffEq.ode113())
"SciPy: RK45" => Dict(:alg=>SciPyDiffEq.RK45())
"SciPy: LSODA" => Dict(:alg=>SciPyDiffEq.LSODA())
"Sundials: Adams" => Dict(:alg=>Sundials.CVODE_Adams())
"liblsoda: LSODA" => Dict(:alg=>LSODA.lsoda())
]
names = first.(_setups)
setups = last.(_setups)
solver_colors = [1 1 1 2 2 3 3 4 4 5 6];
abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10)
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names=names,
print_names=true,
appxsol=test_sol,
dense=DENSE,
save_everystep=SAVE_EVERYSTEP,
numruns=10,
maxiters=Int(1e5),
timeseries_errors=false,
verbose=false
)
Julia: Tsit5 Julia: Vern9 Julia: RadauIIA5 Hairer: dopri5 Hairer: rodas MATLAB: ode45 MATLAB: ode113 SciPy: RK45 SciPy: LSODA Sundials: Adams liblsoda: LSODA
┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead │ caller = npyinitialize() at numpy.jl:67 └ @ PyCall /home/nbosch/.julia/packages/PyCall/L0fLP/src/numpy.jl:67
WorkPrecisionSet of 11 wps
# ProbNumDiffEq.jl solvers:
_setups = [
"PNDE.jl: EK0(4)" => Dict(:alg=>EK0(order=4, smooth=DENSE))
"PNDE.jl: EK1(5)" => Dict(:alg=>EK1(order=5, smooth=DENSE))
]
names = first.(_setups)
setups = last.(_setups)
pn_wp = WorkPrecisionSet(prob,abstols*1e2,reltols*1e2,setups;
names=names,
print_names=true,
appxsol=test_sol,
dense=DENSE,
save_everystep=SAVE_EVERYSTEP,
numruns=10,
maxiters=Int(1e5),
timeseries_errors=false,
verbose=false
)
PNDE.jl: EK0(4) PNDE.jl: EK1(5)
WorkPrecisionSet of 2 wps
wp_plot = plot(wp, title="Lotka-Volterra (non-stiff)", size=(800,500), color=solver_colors, linewidth=8)
plot!(pn_wp, color=:darkorange, markersize=10, linewidth=14)
wp_plot
# Problem definition and reference solution
vdp = @ode_def begin
dy = μ*((1-x^2)*y - x)
dx = 1*y
end μ
prob = ODEProblem(vdp, [0.0, 2.0], (0.0, 6.3), 1e6)
sol = solve(remake(prob, u0=big.(prob.u0), tspan=big.(prob.tspan)),
RadauIIA5(), abstol=1/10^25, reltol=1/10^25)
test_sol = TestSolution(sol)
plot(sol, title="Van-der-Pol Solution", ylim=[-5, 5], legend=false)
# Non-probabilistic solvers
abstols = 1.0 ./ 10.0 .^ (7:14)
reltols = 1.0 ./ 10.0 .^ (4:11)
_setups = [
"Julia: RadauIIA5" => Dict(:alg=>RadauIIA5())
"Julia: Rodas5" => Dict(:alg=>Rodas5())
"Julia: QNDF" => Dict(:alg=>QNDF())
"Hairer: rodas" => Dict(:alg=>ODEInterfaceDiffEq.rodas())
"Hairer: radau" => Dict(:alg=>ODEInterfaceDiffEq.radau())
#"MATLAB: ode23s" => Dict(:alg=>MATLABDiffEq.ode23s()) # Not able to solve VdP!
"MATLAB: ode15s" => Dict(:alg=>MATLABDiffEq.ode15s())
"SciPy: LSODA" => Dict(:alg=>SciPyDiffEq.LSODA())
"SciPy: BDF" => Dict(:alg=>SciPyDiffEq.BDF())
"Sundials: CVODE" => Dict(:alg=>Sundials.CVODE_BDF())
"liblsoda: LSODA" => Dict(:alg=>LSODA.lsoda())
]
names = first.(_setups)
setups = last.(_setups)
solver_colors = [1 1 1 2 2 3 4 4 5 6];
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
names=names,
print_names=true,
appxsol=test_sol,
dense=DENSE,
save_everystep=SAVE_EVERYSTEP,
numruns=10,
maxiters=Int(1e6),
timeseries_errors=false,
verbose=false)
Julia: RadauIIA5 Julia: Rodas5 Julia: QNDF Hairer: rodas Hairer: radau MATLAB: ode15s SciPy: LSODA SciPy: BDF Sundials: CVODE liblsoda: LSODA
WorkPrecisionSet of 10 wps
# ProbNumDiffEq.jl solvers:
_setups = [
"PNDE.jl: EK1(5)" => Dict(:alg=>EK1(order=5, smooth=DENSE))
"PNDE.jl: EK1(7)" => Dict(:alg=>EK1(order=7, smooth=DENSE))
]
names = first.(_setups)
setups = last.(_setups)
pn_wp = WorkPrecisionSet(prob,abstols*1e2,reltols*1e2,setups;
names=names,
print_names=true,
appxsol=test_sol,
dense=DENSE,
save_everystep=SAVE_EVERYSTEP,
numruns=10,
maxiters=Int(1e5),
timeseries_errors=false,
verbose=false
)
PNDE.jl: EK1(5) PNDE.jl: EK1(7)
WorkPrecisionSet of 2 wps
wp_plot = plot(wp, title="Van-der-Pol (stiff)", size=(900,500), color=solver_colors,
linewidth=8, xlims=(1e-12, 1e-1), ylims=(1e-3, 1e1)
)
plot!(pn_wp, color=:darkorange, markersize=10, linewidth=14)
wp_plot
┌ Warning: Invalid negative or zero value 0.0 found at series index 5 for log10 based xscale └ @ Plots /home/nbosch/.julia/packages/Plots/4oFWe/src/utils.jl:95 ┌ Warning: Invalid negative or zero value 0.0 found at series index 6 for log10 based xscale └ @ Plots /home/nbosch/.julia/packages/Plots/4oFWe/src/utils.jl:95